#!/usr/bin/perl -w
use strict;
use DBI;
use DBD::Pg;
use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use uniquify;
use Data::Dump qw(dump);
use db_parameters;
use ortho_species;

print <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_syntenic_distances

		Do a round trip via the neighbour of an ortholog
			and see how far we are from where we started

888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

my $max_cnt_ortholog_sets = 5;
my $max_cnt_gene_pos = 20;

#_________________________________________________________________________________________
#
#	Get the relative gene positions of all orthologs for a species
#
#_________________________________________________________________________________________
sub get_gene_id_vs_pos($$$\%\%)
{
	my ($dbh, $species, $protocol_id, $gene_id_to_pos, $gene_pos_to_id) = @_;

	print STDERR "\t\tGet gene positions for $species:\n";
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
        gene_id, 
        chromosome 
    FROM 
        orthologs.ortholog_sets JOIN 
        taxon.gene_loci USING (gene_id)
    WHERE 
        protocol_id = $protocol_id AND
        species = '$species' AND
        ortholog_id IN 
        (
            SELECT 
                    DISTINCT ON (ortholog_id) 
                    ortholog_id
                FROM
                    orthologs.ortholog_sets
                WHERE 
                    protocol_id = $protocol_id
                GROUP BY 
                    ortholog_id
                HAVING
                    count(DISTINCT species) > 1
        )   
    ORDER BY gene_position;
PL/SQLCMD
	my $curr_chromosome = '';
	my $gene_position = -1;
	for my $data(@{$dbh->selectall_arrayref($sql_cmd)})
	{
		++$gene_position;
		if ($data->[1] ne $curr_chromosome)
		{
			$gene_position += 20000;
			$curr_chromosome = $data->[1];
		}
		$gene_id_to_pos->{$data->[0]} = $gene_position;
		$gene_pos_to_id->{$gene_position} = $data->[0];
	}
	# marker for last gene
	$gene_pos_to_id->{++$gene_position} = undef;

}

#_________________________________________________________________________________________
#
#	get all ortholog sets for a species
#
#_________________________________________________________________________________________
sub get_ortholog_sets($$$\%\%)
{
	my ($dbh, $species, $protocol_id, $gene_id_to_ortholog_id, $ortholog_id_to_gene_ids) = @_;

	print STDERR  "\t\tGet ortholog sets for $species:\n";
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
        gene_id, 
        ortholog_id
    FROM 
        orthologs.ortholog_sets
    WHERE 
        protocol_id = $protocol_id  AND
        species = '$species'  AND
        ortholog_id IN 
        (
            SELECT 
                    DISTINCT ON (ortholog_id) 
                    ortholog_id
                FROM
                    orthologs.ortholog_sets
                WHERE 
                    protocol_id = $protocol_id
                GROUP BY 
                    ortholog_id
                HAVING
                    count(DISTINCT species) > 1
        )
PL/SQLCMD
# and        ortholog_type ~ 'to'
	for my $data(@{$dbh->selectall_arrayref($sql_cmd)})
	{
		my ($gene_id,  $ortholog_id) = @$data;
		$gene_id_to_ortholog_id->{$gene_id} = $ortholog_id;
		push(@{$ortholog_id_to_gene_ids->{$ortholog_id}},  $gene_id);
	}
}


my %promiscuous_ortholog_ids;
#_________________________________________________________________________________________
#
#	get_ortholog_ids_for_neighbours_of_ortholog
#
#
#		for each ortholog in the other species, find neighbours
#			for each neighbour find back orthologs in the original species
#		
#_________________________________________________________________________________________
sub get_ortholog_ids_for_neighbours_of_ortholog($$\%\%\%$)
{
	my ($ortholog_other_gene_id, $direction, $gene_id_to_pos, 
		$gene_id_to_ortholog_id, $other_gene_pos_to_id, $max_cnt_ortholog_sets) =@_;

	# get gene position of neighbour of current ortholog
	die "$ortholog_other_gene_id" unless exists $gene_id_to_pos->{$ortholog_other_gene_id};
	my $other_curr_gene_pos	= $gene_id_to_pos->{$ortholog_other_gene_id};
	my $ortholog_id			= $gene_id_to_ortholog_id->{$ortholog_other_gene_id};

	# ortholog ids of neighbours. Collect max_cnt of these
	my %ortholog_ids;
	++$ortholog_ids{$ortholog_id};

	#DEBUG
	#my $starting_pos = $other_curr_gene_pos;
	
	my $gene_pos_dist = 0;
	while (1)
	{
		++$gene_pos_dist;
		# stop if already have enough ortholog ids
		if (keys %ortholog_ids > $max_cnt_ortholog_sets && $gene_pos_dist > $max_cnt_gene_pos)
		{
			last;
		}

		#
		# stop on interchromosome gap
		$other_curr_gene_pos += $direction;
		unless (exists $other_gene_pos_to_id->{$other_curr_gene_pos} &&
				defined $other_gene_pos_to_id->{$other_curr_gene_pos})
		{
			#print STDERR "\n\n",dump(%ortholog_ids ), "\n";							#DEBUG
			#print STDERR "starting_pos = $starting_pos\n";								#DEBUG
			#print STDERR $gene_id_to_pos->{$ortholog_other_gene_id}, "\n";				#DEBUG
			#print STDERR "$ortholog_other_gene_id $direction: $other_curr_gene_pos";	#DEBUG
			last;
		}

		
		my $gene_id_of_current_neighbour = $other_gene_pos_to_id->{$other_curr_gene_pos};
		my $neighbour_ortholog_id = $gene_id_to_ortholog_id->{$gene_id_of_current_neighbour};

		# ignore this ortholog if this is part of a 'promiscuous' orthologues set
		# with links all over the genome
		next if (exists $promiscuous_ortholog_ids{$neighbour_ortholog_id});

		# save this ortholog_id
		++$ortholog_ids{$neighbour_ortholog_id};
	}
	
	# remove starting ortholog
	delete ($ortholog_ids{$ortholog_id});
	#print STDERR "!";	#DEBUG
	return keys %ortholog_ids;	
}


#_________________________________________________________________________________________
#
#	min
#		
#_________________________________________________________________________________________
sub min
{
	return () unless @_;
	my @array = sort {$a <=> $b} (@_);
	return $array[0];
}

#_________________________________________________________________________________________
#
#	get_minimum_dist
#
#		from one gene to any in list
#		
#_________________________________________________________________________________________
sub get_minimum_dist(\@$\%)
{
	my ($neighbour_gene_ids, $gene_id, $gene_id_to_pos) = @_;
	return 100000000 unless (@$neighbour_gene_ids);
	# for debug
	#for (@$neighbour_gene_ids)
	#{
	#	die unless $_;
	#	die "$_" unless exists $gene_id_to_pos->{$_};
	#}
	for (@$neighbour_gene_ids) {die $_ unless exists $gene_id_to_pos->{$_} } 
	my @neighbour_gene_pos = map {$gene_id_to_pos->{$_}} @$neighbour_gene_ids;
	my $gene_pos = $gene_id_to_pos->{$gene_id};
	die "$gene_id" unless defined $gene_pos;
	return min(20000, map{abs($_ - $gene_pos)} @neighbour_gene_pos);
}


my %ignore_gene_ids;
sub gene_ids_from_ortholog_ids(\@\%)
{
	my ($ortholog_ids, $ortholog_id_to_gene_ids) = @_;
	my @gene_ids;
	for my $ortholog_id(@$ortholog_ids)
	{
		for (@{$ortholog_id_to_gene_ids->{$ortholog_id}})
		{
			push(@gene_ids, $_) unless exists $ignore_gene_ids{$_};
		}
	}
	return @gene_ids;
}
sub get_syntenic_distance_in_single_direction(\%\%\@$$\%\%\%)
{
	my ($my_ortholog_id_to_gene_ids,
		$other_gene_pos_to_id,
		$ortholog_other_gene_ids, 
		$direction,
		$ortholog_id,
		$gene_id_to_pos, 
		$gene_id_to_ortholog_id, 
		$gene_id_to_dist) = @_;
	# get all up stream neighbours of my orthologs in the other species
	my @back_orthologs_ids;

	for my $ortholog_other_gene_id(@$ortholog_other_gene_ids)
	{
		push(@back_orthologs_ids, 
				get_ortholog_ids_for_neighbours_of_ortholog($ortholog_other_gene_id, 
															$direction, 
															%$gene_id_to_pos, 
															%$gene_id_to_ortholog_id, 
															%$other_gene_pos_to_id,
															$max_cnt_ortholog_sets));
	}

	#
	#	convert ortholog_ids into gene ids for my own species
	#
	@back_orthologs_ids = uniquify(@back_orthologs_ids);
	if (!@back_orthologs_ids)
	{
		for my $gene_id(@{$my_ortholog_id_to_gene_ids->{$ortholog_id}})
		{
			$gene_id_to_dist->{$gene_id}[($direction > 0) ? 0: 1] = 20000;
		}
		return;
	}

	my @back_orthologs_gene_ids = gene_ids_from_ortholog_ids(@back_orthologs_ids,
															 %$my_ortholog_id_to_gene_ids);
	return unless @back_orthologs_gene_ids;

	#
	#	Get minimum syntenic distances per ortholog set and per gene
	#
	for my $gene_id(@{$my_ortholog_id_to_gene_ids->{$ortholog_id}})
	{
		# get syntenic distances
		my $dist = get_minimum_dist( @back_orthologs_gene_ids, 
									 $gene_id,
									 %$gene_id_to_pos);
		$gene_id_to_dist->{$gene_id}[($direction > 0) ? 0: 1] = $dist;
	}
}

#_________________________________________________________________________________________
#
#	get_promiscuous_orthologs
#
#		Get ortholog sets whcih have greater than 4 chromosomes in either species
#		We do a orthology round trip via the neighbours in the other species and
#			see how far away we are from where we started
#		
#_________________________________________________________________________________________
sub get_promiscuous_orthologs($\%)
{
	my ($dbh, $promiscuous_ortholog_ids) = @_;

	print STDERR "\t\tGet promiscuous orthologs:\n";
	my $sql_cmd = <<"PL/SQLCMD";
	SELECT 
			DISTINCT ON (ortholog_id) 
			ortholog_id
		FROM 
			orthologs.ortholog_sets NATURAL JOIN 
			taxon.gene_loci 
		WHERE 
			protocol_id = $curr_protocol_id
		GROUP BY 
			ortholog_id, 
			species 
		HAVING 
			COUNT(chromosome) > 4;
PL/SQLCMD
	for my $data(@{$dbh->selectall_arrayref($sql_cmd)})
	{
		++$promiscuous_ortholog_ids->{$data->[0]};
	}
}


#_________________________________________________________________________________________
#
#	get_syntenic_distances
#
#		for each gene and each ortholog set find the minimum syntenic distance
#		We do a orthology round trip via the neighbours in the other species and
#			see how far away we are from where we started
#		1) go to ortholog in other species
#		2) get neighbouring X number of ortholog sets (distinct ortholog_ids)
#		3) Find the orthologs for those in my own species
#		
#_________________________________________________________________________________________
sub get_syntenic_distances(\%\%\%\%\%\%)
{
	my ($my_ortholog_id_to_gene_ids, 		
		$other_ortholog_id_to_gene_ids, 
		$other_gene_pos_to_id, 
		$gene_id_to_pos, 
		$gene_id_to_ortholog_id,
		$gene_id_to_dist) = @_;
	print STDERR "\t", scalar keys %$my_ortholog_id_to_gene_ids, " orthologue sets\n\t[";
	my $cnt_ortholog_sets = 0;
	my $every = (scalar keys %$my_ortholog_id_to_gene_ids) / 200 + 1;

	for my $ortholog_id(keys %$my_ortholog_id_to_gene_ids)
	{
		if (++$cnt_ortholog_sets % $every == 0)
		{
			my $line = int($cnt_ortholog_sets / $every);
			print STDERR "\n\t" if ($line && $line % 50 == 0);
			print STDERR ".";
		}
		# get all orthologs
		my @ortholog_other_gene_ids = @{$other_ortholog_id_to_gene_ids->{$ortholog_id}};
		
		get_syntenic_distance_in_single_direction(%$my_ortholog_id_to_gene_ids,
												  %$other_gene_pos_to_id,
												  @ortholog_other_gene_ids, 
												  1,
												  $ortholog_id,
												  %$gene_id_to_pos, 
												  %$gene_id_to_ortholog_id, 
												  %$gene_id_to_dist);

		get_syntenic_distance_in_single_direction(%$my_ortholog_id_to_gene_ids,
												  %$other_gene_pos_to_id,
												  @ortholog_other_gene_ids, 
												  -1,
												  $ortholog_id,
												  %$gene_id_to_pos, 
												  %$gene_id_to_ortholog_id, 
												  %$gene_id_to_dist);
		
	}
	print STDERR "]\n";
}

#
#	Connect to genome database
#
my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);

#
#	get_promiscuous_orthologs
#
#		use global: yuck
get_promiscuous_orthologs($dbh, %promiscuous_ortholog_ids);




#DEBUG_1
$dbh->do("DELETE FROM orthologs.syntenic_distances where protocol_id = $curr_protocol_id;");
$dbh->do("VACUUM ANALYSE orthologs.syntenic_distances;");
for my $protocol (get_orthologs_protocols())
{
	my ($protocol_name,  $protocol_id) = @$protocol;
	print STDERR "\tOrthology protocol [$protocol_name]...\n";
	

	#DEBUG_1
	#next unless $protocol_id == 1;
	next unless $protocol_id == $curr_protocol_id;

	my (%gene_id_to_pos, %gene_pos_to_id1, %gene_pos_to_id2);
	my (%gene_id_to_ortholog_id, %ortholog_id_to_gene_ids1, %ortholog_id_to_gene_ids2);

	#
	#	get gene positions to and from gene_id
	#
	get_gene_id_vs_pos($dbh, $ortho_name1, $protocol_id, 
											%gene_id_to_pos, %gene_pos_to_id1);
	get_gene_id_vs_pos($dbh, $ortho_name2, $protocol_id, 
											%gene_id_to_pos, %gene_pos_to_id2);

	#
	#	get ortholog_id to and from gene_id
	#
	get_ortholog_sets($dbh, $ortho_name1, $protocol_id, %gene_id_to_ortholog_id, 
											%ortholog_id_to_gene_ids1);
	get_ortholog_sets($dbh, $ortho_name2, $protocol_id, %gene_id_to_ortholog_id, 
											%ortholog_id_to_gene_ids2);


	#
	#	For each gene find distance to its closest synteny partner via its ortholog in
	#		the other species
	#
	my (%gene_id_to_dist);
	print STDERR "\t\tCalculating orthologue syntenic distances starting from $ortho_name1...\n";
	get_syntenic_distances(	%ortholog_id_to_gene_ids1, %ortholog_id_to_gene_ids2, %gene_pos_to_id2, 
							%gene_id_to_pos, %gene_id_to_ortholog_id,
							%gene_id_to_dist);
	print STDERR "\t\tCalculating orthologue syntenic distances starting from $ortho_name2...\n";
	get_syntenic_distances(	%ortholog_id_to_gene_ids2, %ortholog_id_to_gene_ids1, %gene_pos_to_id1, 
							%gene_id_to_pos, %gene_id_to_ortholog_id,
							%gene_id_to_dist);

	#
	#	run get syntenic distances again with the non-syntenic genes for the first run on
	#		a black list
	#
	while( my ($gene_id,  $dist) = each %gene_id_to_dist)
	{
		my @dist = sort {$a <=> $b } (@$dist);
		$ignore_gene_ids{$gene_id}++ if ($dist[0] > 20);
	}

	print STDERR "\t\tCalculating orthologue syntenic distances starting from $ortho_name1...\n";
	get_syntenic_distances(	%ortholog_id_to_gene_ids1, %ortholog_id_to_gene_ids2, %gene_pos_to_id2, 
							%gene_id_to_pos, %gene_id_to_ortholog_id,
							%gene_id_to_dist);
	print STDERR "\t\tCalculating orthologue syntenic distances starting from $ortho_name2...\n";
	get_syntenic_distances(	%ortholog_id_to_gene_ids2, %ortholog_id_to_gene_ids1, %gene_pos_to_id1, 
							%gene_id_to_pos, %gene_id_to_ortholog_id,
							%gene_id_to_dist);


	print STDERR "\t\tEntering into panda...\n";
	$dbh->do("COPY orthologs.syntenic_distances(".
					"gene_id, closest_neighbour_dist, ".
					"other_dir_neighbour_dist, protocol_id ) FROM STDIN") or die;
	my $cnt_entries = 0;
	while( my ($gene_id,  $dist) = each %gene_id_to_dist)
	{
		$dbh->pg_putline(join ("\t",  $gene_id,  (sort {$a <=> $b }(@$dist)), $protocol_id). "\n");
	}
	$dbh->pg_endcopy;

	#clear black list
	%ignore_gene_ids = ();
}


	#_____________________________________________________________________________________
	#
	#	If filtering on orthologs, add syntenic distances for pseudos
	#		
	#_____________________________________________________________________________________
	print STDERR "\tSet unknown chromosomes to have a syntenic distance of 100000000...\n";
		$dbh->do( <<"PL/SQLCMD");
        INSERT INTO 
            orthologs.syntenic_distances
            (gene_id, 
             protocol_id, 
             closest_neighbour_dist, 
             other_dir_neighbour_dist) 
            SELECT 
                    gene_id, 
                    protocol_id, 
                    20000,
                    20000
                FROM
                    orthologs.syntenic_distances sd RIGHT JOIN 
                    orthologs.ortholog_sets USING (gene_id, protocol_id)
                WHERE
                    sd.gene_id IS NULL;
PL/SQLCMD
	#_____________________________________________________________________________________
	#
	#	Set unknown chromosomes to have a syntenic distance of 100000000...
	#		
	#_____________________________________________________________________________________
	print STDERR "\tSet unknown chromosomes to have a syntenic distance of 100000000...\n";
	{
		$dbh->do( <<"PL/SQLCMD");
        UPDATE 
				orthologs.syntenic_distances 
			SET 
				closest_neighbour_dist = 50000 
            FROM 
                taxon.gene_loci
			WHERE 
				closest_neighbour_dist >= 20000 AND 
				gene_loci.gene_id = syntenic_distances.gene_id AND NOT 
				gene_loci.chromosome ~ '^[A-Z0-9]+\$';
		UPDATE 
				orthologs.syntenic_distances 
			SET 
				other_dir_neighbour_dist = 50000 
            FROM 
                taxon.gene_loci
			WHERE 
				other_dir_neighbour_dist >= 20000 AND 
				gene_loci.gene_id = syntenic_distances.gene_id AND NOT 
				gene_loci.chromosome ~ '^[A-Z0-9]+\$';
PL/SQLCMD
	}
	print STDERR "\tDone\n";
	print STDERR "\tCompleted\n";
	log_pipeline_stage($dbh, $curr_protocol_id, 34, 'table syntenic_distances', $start_time);
	$dbh->disconnect();

